Spatial Statistics 1 - Notes

UK administrative Geographies

The UK consists of four countries: England, Scotland, Wales and Norther Ireland (in order of population size). Great Britain refers to the island that has the first three countries. Administrative data for the UK often follows country divides, although England and Wales are accounted for together for many things by the Office of National Statistics. Administrative geographies are one of 8 overarching types of statistical geographies, as seen below:

Broadly, administrative geographies are most general purpose for the UK as this is the level at which policy is made. Also, the units of the Statistical Building Blocks constrained are by these. For example, middle-level output areas (MSOAs) and the smaller geographies below all aggregate to Local Authority Districts. Administrative geographies come in several levels of detail:

Worth noting about the Statistical Building Blocks is that they are derived from populations counts, not areas. Below is an overview of the thresholds used to create these geographies.

More about these population-weighted geographies here

Getting Data

Polygons

You can get geographic data for the UK from the open geography portal via an API call. This is convenient because it means you don’t have to store large files on your machine and can share your work easier. For simplicity we will use regions for example below. There are nine regions in England.

library(raster)
library(knitr)
library(geojsonio)

library(sp)
library(tmap)

library(spdep)

library(reshape2)

library(rsq)

#connect to the open geography portal API
regions_json <- geojson_read("https://opendata.arcgis.com/datasets/8d3a9e6e7bd445e2bdcc26cdf007eac7_1.geojson", what = "sp")
z_regions_json <- regions_json
plot(regions_json )

Variables

I’m using table WU02EW - Location of usual residence and place of work by age. This data is available down to MSOA, but we will be using it at regional level for England.

#Connect to the NOMIS API to get Data
#Note: I'm only using England here for convience.
t_wu02Ew  <- read.csv(file = "https://www.nomisweb.co.uk/api/v01/dataset/NM_1206_1.data.csv?date=latest&usual_residence=2013265921...2013265930&place_of_work=2013265921...2013265930&age=0...6&measures=20100", header=TRUE)
kable(head(t_wu02Ew))

DATE	DATE_NAME	DATE_CODE	DATE_TYPE	DATE_TYPECODE	DATE_SORTORDER	USUAL_RESIDENCE	USUAL_RESIDENCE_NAME	USUAL_RESIDENCE_CODE	USUAL_RESIDENCE_TYPE	USUAL_RESIDENCE_TYPECODE	USUAL_RESIDENCE_SORTORDER	PLACE_OF_WORK	PLACE_OF_WORK_NAME	PLACE_OF_WORK_CODE	PLACE_OF_WORK_TYPE	PLACE_OF_WORK_TYPECODE	PLACE_OF_WORK_SORTORDER	AGE	AGE_NAME	AGE_CODE	AGE_TYPE	AGE_TYPECODE	AGE_SORTORDER	MEASURES	MEASURES_NAME	OBS_VALUE	OBS_STATUS	OBS_STATUS_NAME	OBS_CONF	OBS_CONF_NAME	URN	RECORD_OFFSET	RECORD_COUNT
2011	2011	2011	date	0	0	2013265921	North East	E12000001	regions	480	0	2013265921	North East	E12000001	regions	480	0	0	Aged 16 and over	0	Age	1000	0	20100	Value	936525	A	Normal Value	FALSE	Free (free for publication)	Nm-1206d1d32176e1d2013265921d2013265921d0d20100	0	700
2011	2011	2011	date	0	0	2013265921	North East	E12000001	regions	480	0	2013265921	North East	E12000001	regions	480	0	1	Aged 16-24	1	Age	1000	1	20100	Value	130827	A	Normal Value	FALSE	Free (free for publication)	Nm-1206d1d32176e1d2013265921d2013265921d1d20100	1	700
2011	2011	2011	date	0	0	2013265921	North East	E12000001	regions	480	0	2013265921	North East	E12000001	regions	480	0	2	Aged 25-34	2	Age	1000	2	20100	Value	199584	A	Normal Value	FALSE	Free (free for publication)	Nm-1206d1d32176e1d2013265921d2013265921d2d20100	2	700
2011	2011	2011	date	0	0	2013265921	North East	E12000001	regions	480	0	2013265921	North East	E12000001	regions	480	0	3	Aged 35-49	3	Age	1000	3	20100	Value	344176	A	Normal Value	FALSE	Free (free for publication)	Nm-1206d1d32176e1d2013265921d2013265921d3d20100	3	700
2011	2011	2011	date	0	0	2013265921	North East	E12000001	regions	480	0	2013265921	North East	E12000001	regions	480	0	4	Aged 50-64	4	Age	1000	4	20100	Value	243426	A	Normal Value	FALSE	Free (free for publication)	Nm-1206d1d32176e1d2013265921d2013265921d4d20100	4	700
2011	2011	2011	date	0	0	2013265921	North East	E12000001	regions	480	0	2013265921	North East	E12000001	regions	480	0	5	Aged 65-74	5	Age	1000	5	20100	Value	15701	A	Normal Value	FALSE	Free (free for publication)	Nm-1206d1d32176e1d2013265921d2013265921d5d20100	5	700

This dataset is currently a flat 2D table despite containing multi-way tabulated counts: by region,origin,destination and age group. To wrangle this into counts by geographic unit to use with the polygon data, we apply the following transformation:

#names(t_wu02Ew )#Select only columns we need for now
ac <- c("USUAL_RESIDENCE_CODE","AGE_NAME","OBS_VALUE")
t_wu02Ew_ac <-t_wu02Ew[,ac]

reg_var <- xtabs(OBS_VALUE ~ USUAL_RESIDENCE_CODE + AGE_NAME, data=t_wu02Ew_ac)
reg_var <- as.data.frame.matrix(reg_var) 
kable(reg_var)

	Aged 16-24	Aged 16 and over	Aged 25-34	Aged 35-49	Aged 50-64	Aged 65-74	Aged 75+
E12000001	137047	974625	208221	358197	252002	16252	2906
E12000002	386853	2705931	598669	986277	666548	57777	9807
E12000003	294930	2029907	442170	739320	505157	41209	7121
E12000004	247828	1777612	372762	656200	454732	39633	6457
E12000005	290056	2106075	458558	771101	526138	51665	8557
E12000006	315871	2299955	496908	834637	581497	60791	10251
E12000007	374191	3197606	1056357	1101460	591100	60959	13539
E12000008	461335	3391170	731448	1238930	852032	91961	15464
E12000009	292428	2024395	417547	716798	531510	56730	9382
W92000004	160618	1117784	238664	404170	284858	25071	4403

This data gives information on the working population by place of residence by age group. For simplicity, we will combine the seven age groups into two. Let’s assume person over 65 could be retired, whilst all other ages we can expect to working to generate two variables:

retired <- c("Aged 65-74", "Aged 75+" )
reg_var$w_Age <- rowSums(reg_var[,!names(reg_var) %in% retired])
reg_var$r_Age <- rowSums(reg_var[,retired])

Combining the two

Although you could join your tables on region names, geography codes where available will give you a much cleaner merge.

#attach the data to the dataframe component of the Spatial data
#Note: 0 means attach by row names
regions_json@data <- merge(regions_json@data, reg_var, by.x= "rgn15cd", by.y=0 )

Plotting Your Maps

A useful library for showing your geographic data is tmap. Dr. Robin Lovelace has a good primer on using this as you can see in the link above. Notice from the map below that we have more observations of individuals working above the age of 65 in the south and south-west.

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(regions_json) + tm_polygons(col="r_Age")

So far we have covered the basics of UK geographies and data sources, including how to call and wrangle your data.

---

Spatial Autocorrelation

This happens when observations that are close to each other show some dependence - i.e. close observations can to some extent predict each other and show a spatial patterning. This self-correlation can be a problem because it violates the assumption of independence some methods rely one. The main two point in flagging up such a spatial dependence are: firstly, defining who are your neighbours; and secondly, what your relationship to these is (which are expressed as weights).

One common way to define neighbours is by contiguity. This means sharing a common boundary and you can generally check this by checking for at least one common boundary point. This neighbouring definition is what I will use for the examples below. However note that neighbours can equally be defined by some more sophisticated network measure, as well a social or physical similarity between places.

Moran’s I

Theory

The most common way to measure spatial autocorrelation is using Moran’s I, specificially Global Moran’s I.
\[ I=\frac{N}{W} \frac{\sum_{i} \sum_{j} w_{i j}\left(x_{i}-\overline{x}\right)\left(x_{j}-\overline{x}\right)}{\sum_{i}\left(x_{i}-\overline{x}\right)^{2}} \] where \({\displaystyle N}\) is the number of spatial units indexed by \({\displaystyle i}\) i and \({\displaystyle j}\) j
\({\displaystyle x}\) is the variable of interest
\({\displaystyle {\bar {x}}}\) is the mean of \({\displaystyle x}\)
\({\displaystyle w_{ij}}\) is a matrix of spatial weights with zeroes on the diagonal (i.e., \({\displaystyle w_{ii}=0}\) \({\displaystyle w_{ii}=0}\))
\({\displaystyle W}\) is the sum of all \({\displaystyle w_{ij}}\) \(w_{ij}\)

This is an index that ranges from -1 to 1 indicating whether observations are disperse or clustered. To decide if the variance of your variable in space is actually significant, Moran’s I is often reported with a p-value assuming the null-hypothesis that observations are independent in space.

Disperions images source.

Example

Defining your neighbours

#queen=T; neighbour can share only point
neighbs<- poly2nb(regions_json, queen=TRUE)
summary(neighbs)

## Neighbour list object:
## Number of regions: 9 
## Number of nonzero links: 30 
## Percentage nonzero weights: 37.03704 
## Average number of links: 3.333333 
## Link number distribution:
## 
## 2 3 4 5 
## 3 2 2 2 
## 3 least connected regions:
## 0 6 8 with 2 links
## 2 most connected regions:
## 3 7 with 5 links

Notice that we have links created for all nine regions. These can be access by:

regions_json$rgn15nm[1] # region

## [1] North East
## 9 Levels: East Midlands East of England London North East ... Yorkshire and The Humber

regions_json$rgn15nm[unlist(neighbs[[1]])] # its neighbours

## [1] North West               Yorkshire and The Humber
## 9 Levels: East Midlands East of England London North East ... Yorkshire and The Humber

Establishing your relationship to these
Now that the neighbours are defined, their relation needs to be articulated. This is done with the style param. The simplest way, which I will use here, is to assume equal weight for all neighbours:

lw <- nb2listw(neighbs, style="W")

Moran’s Global I
Having defined the neighbours and their relations, we can check our index and its corresponding p-value in two ways:

moran.test(regions_json@data$r_Age,lw)

## 
##  Moran I test under randomisation
## 
## data:  regions_json@data$r_Age  
## weights: lw    
## 
## Moran I statistic standard deviate = 1.5783, p-value = 0.05725
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.16636716       -0.12500000        0.03407958

#Permutation test for Moran's I statistic
MC<- moran.mc(regions_json@data$r_Age, lw, nsim=599)
MC

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  regions_json@data$r_Age 
## weights: lw  
## number of simulations + 1: 600 
## 
## statistic = 0.16637, observed rank = 544, p-value = 0.09333
## alternative hypothesis: greater

Note the above provides same index but different p-values. In both cases, the spatial distribution of our retired group is slightly correlated in space. Specifically this group is slightly clustered with a significance level of p > 0.05.

The different p-values arise because the two methods make different assumptions about the null hypothesis, i.e. what spatially independent actually means. At a higher level, it has to do with what assumption we make about the variance of our variable. In the first example we are assuming that no spatial dependency means that each value is likely to occur in any polygon. We do something similar in the second example, but instead of looking at all polygons at once, we make this assumption iteratively on subsets of our data to define our distribution. Below is a plot of what the simulated distribution looks like with the line indicating the global value of \(I\).

plot(MC, main="", las=1)

When dealing with many polygons a simulation-based approach would make more sense. Brunsdon & Comber go into this with more detail.. There are additional ways to visualise spatial autocorrelation, notably using variograms. You can find a useful starting point with code examples here.)

Further Reading on Moran’s I

Useful overview from ArcGIS Manuel Gimond’s extension to local Morans’ I Manual for spdep

Geary’s C

A less popular alternative to Moran’s I is Geary’s C, also known as Geary’s contiguity ratio or simply Geary’s ratio. It ranges from 0 to 2. Inversely to \(I\), lower values \(G\) indicate clustering and higher values indicate dispersal. I looks very similar to \(I\):

\[ C=\frac{(N-1) \sum_{i} \sum_{j} w_{i j}\left(x_{i}-x_{j}\right)^{2}}{2 W \sum_{i}\left(x_{i}-\overline{x}\right)^{2}} \]

where \({\displaystyle N}\) is the number of spatial units indexed by \({\displaystyle i}\) i and \({\displaystyle j}\) j
\({\displaystyle x}\) is the variable of interest
\({\displaystyle {\bar {x}}}\) is the mean of \({\displaystyle x}\)
\({\displaystyle w_{ij}}\) is a matrix of spatial weights with zeroes on the diagonal (i.e., \({\displaystyle w_{ii}=0}\) \({\displaystyle w_{ii}=0}\))
\({\displaystyle W}\) is the sum of all \({\displaystyle w_{ij}}\) \(w_{ij}\)

In the case of Moran’s \(I\), the numerator considers the value per polygon for each neighbour relative to the mean \({\displaystyle {\bar {x}}}\). By contrast, in Geary’s \(C\), we consider the relationship between spatial units directly with \(\left(x_{i}-x_{j}\right)^{2}\), penalising large difference by squaring the whole expression. As we are working with two units at a time, how we normalise also changes: we do this using \({\displaystyle N}\) - 1 units and double the sum of all weights. As \(C\) relies less on the mean of all observations which shifts drastically with extreme outliers, it provides a more nuanced perspective on local variations.

geary.test(regions_json@data$r_Age,lw)

## 
##  Geary C test under randomisation
## 
## data:  regions_json@data$r_Age 
## weights: lw 
## 
## Geary C statistic standard deviate = 0.87315, p-value = 0.1913
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.82089773        1.00000000        0.04207492

#Permutation test for Moran's I statistic
MC<- geary.mc(regions_json@data$r_Age, lw, nsim=599)
MC

## 
##  Monte-Carlo simulation of Geary C
## 
## data:  regions_json@data$r_Age 
## weights: lw 
## number of simulations + 1: 600 
## 
## statistic = 0.8209, observed rank = 111, p-value = 0.185
## alternative hypothesis: greater

Further Examples using Geary’s C

Getis and Ord’s G

At this point I have considered two global measures. These can be extended and calculated on a smaller subset to produce local measures. However, an intrinsically local measure is Getis and Ord’s \(G\), which gives an indication of clusters based on their values:

\[ G=\frac{\sum_{i=1} \sum_{j=1} w_{i, j} x_{i} x_{j}}{\sum_{i=1} \sum_{j=1} x_{i} x_{j}} \]

where the number spatial units indexed by \({\displaystyle i}\) i and \({\displaystyle j}\) j
\({\displaystyle x}\) is the variable of interest
\({\displaystyle w_{ij}}\) is a matrix of spatial weights with zeroes on the diagonal (i.e., \({\displaystyle w_{ii}=0}\) \({\displaystyle w_{ii}=0}\))

Notice, that in the above units are considered relative to each other and their relationship (ie. weights) only. Getis and Ord’s \(G\) is useful for identifying clusters based on their values, low or high, and is therefore good for identifying hotspots.

You can see worked \(G\) example here.

Finally, for consistency, I will just mention a final measure Ripleys’s K. This is a descriptive spatial statistic that shows autocorrelation based on multiple distance bands. It is different from all the above because it is not being strictly bound to the whole system or just a units’ local neighbours, and so is useful to see how spatial patterns change with scale. Dr. David Murrel will go over this shortly.

Thus, we have an overview of the core descriptive spatial statistics.

---

Gravity Models

Gravity models are used for estimating movements between two places. These have be used in the context of migration and in transport for trip estimations. Given the wide range of applications, you can guess that this is a big topic with many variations on it. What all of these have in common, is that they are based on Newton’s law of gravitation:

\[ F_{i j}=G \frac{M_{i} M_{j}}{D_{i j}^{2}} \] where \(M_{i}\) and \(M_{j}\) are the mass of two objects \({\displaystyle G}\) is a constant that denotes gravitational pull
\({D_{i j}}\) is the distance between the two objects

For movement modelling, Newton’s force is adapted to describe the flow between two places. Masses are replaced by area characteristics: where you have some constant multiplied by some observations at the origins \(V_{i}^{\mu}\) and at the destinations \(W_{j}^{\alpha}\). In the numerator you have distance, which basically says that the force diminishes over space.

\[ T_{i j}= \frac{V_{i}^{\mu} W_{j}^{\omega}}{d_{i j}^{\beta}} \]

\({D_{i j}}\) generally looks like a power: \[ f\left(d_{i j}\right)=d_{i j}^{\mathcal{B}} \] Or exponential function:

\[ f\left(d_{i j}\right)=\exp \left(\beta d_{i j}\right) \] The shape and choice of the constant \(\beta\) will depend on how much value you give to proximity in predicting flows. Assuming we choose an exponential function, the unconstrainted Gravity equations can be expressed as:

\[ T_{i j}=V_{i}^{\mu} W_{j}^{\alpha} f\left(d_{i j}\right) \] or more conveniently (assuming an exponential distance function) as: \[ \ln T_{i j}=k+\mu \ln V_{i}+\alpha \ln W_{j}-\beta d_{i j} \] Which can be modelled as a log-linear regression model. The choice of origin variables which would make it ‘emissive’ or destination variables with would make places ‘attractive’ is context dependent. However, what all Gravity Models have in common is the need to describe distance between two places (most commonly Euclidean).

Getting Distance

crs(regions_json)

## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

Firstly to get the distance, we need to check our coordinate system. A common mistake is to happily calculate the distance between polygons without verifying units. It is a dangerous one because in relative terms, the distances might look alright, but they won’t correspond to practical world measures. See a baseline for the UK below:

Conveniently, for the BNG projection in the UK, Euclidean distance maps to meters.

#this is the BNG: CRS("+init=epsg:27700")
regions_json <-spTransform(regions_json, CRS("+init=epsg:27700"))
dist_m <- spDists(regions_json)
colnames(dist_m) <- regions_json$rgn15cd
rownames(dist_m) <- regions_json$rgn15cd
kable(dist_m)

	E12000001	E12000002	E12000003	E12000004	E12000005	E12000006	E12000007	E12000008	E12000009
E12000001	0.0	121184.53	125638.67	243896.0	283721.2	348467.20	409329.80	425126.59	454548.3
E12000002	121184.5	0.00	99867.48	179046.7	176972.4	297039.67	333943.03	341419.79	339525.9
E12000003	125638.7	99867.48	0.00	118641.7	179350.2	224625.87	283966.05	301059.44	354195.8
E12000004	243896.0	179046.68	118641.72	0.0	111376.3	118317.43	165433.86	182984.30	267949.6
E12000005	283721.2	176972.43	179350.17	111376.3	0.0	193444.75	184072.14	179457.38	174904.7
E12000006	348467.2	297039.67	224625.87	118317.4	193444.8	0.00	94220.43	128473.25	290312.9
E12000007	409329.8	333943.03	283966.05	165433.9	184072.1	94220.43	0.00	36136.67	219008.5
E12000008	425126.6	341419.79	301059.44	182984.3	179457.4	128473.25	36136.67	0.00	187368.1
E12000009	454548.3	339525.89	354195.80	267949.6	174904.7	290312.89	219008.52	187368.14	0.0

For further use, you can convert the distances into list form.

#tranform the nxn matrix into a pairs list with nxn rows:
dist_pl <- melt(dist_m)
names(dist_pl) <- c('Origin','Destination','Distance_m2')
kable(head(dist_pl))

Origin	Destination	Distance_m2
E12000001	E12000001	0.0
E12000002	E12000001	121184.5
E12000003	E12000001	125638.7
E12000004	E12000001	243896.0
E12000005	E12000001	283721.2
E12000006	E12000001	348467.2

With the theoretical grounding and distance calculation example from above you have all you need to explore different families of gravity models. The most basic is the unconstrained (log-linear Poisson regression model) version as mentioned above.

od <- c("USUAL_RESIDENCE_NAME","USUAL_RESIDENCE_CODE", "PLACE_OF_WORK_NAME","PLACE_OF_WORK_CODE", "AGE_NAME","OBS_VALUE")
t_wu02Ew_od <-t_wu02Ew[,od]
t_wu02Ew_od <-t_wu02Ew_od[t_wu02Ew$AGE_NAME %in% retired,] #get those "Aged 65-74", "Aged 75+"

t_wu02Ew_od <- aggregate(. ~USUAL_RESIDENCE_CODE + PLACE_OF_WORK_CODE + USUAL_RESIDENCE_NAME + PLACE_OF_WORK_NAME, data=t_wu02Ew_od, sum, na.rm=TRUE)
t_wu02Ew_od$AGE_NAME <- NULL

kable(t_wu02Ew_od[1:20,])

USUAL_RESIDENCE_CODE	PLACE_OF_WORK_CODE	USUAL_RESIDENCE_NAME	PLACE_OF_WORK_NAME	OBS_VALUE
E12000006	E12000006	East	East	62825
E12000004	E12000006	East Midlands	East	788
E12000007	E12000006	London	East	2576
E12000001	E12000006	North East	East	34
E12000002	E12000006	North West	East	73
E12000008	E12000006	South East	East	1039
E12000009	E12000006	South West	East	127
W92000004	E12000006	Wales	East	43
E12000005	E12000006	West Midlands	East	94
E12000003	E12000006	Yorkshire and The Humber	East	81
E12000006	E12000004	East	East Midlands	430
E12000004	E12000004	East Midlands	East Midlands	41618
E12000007	E12000004	London	East Midlands	111
E12000001	E12000004	North East	East Midlands	44
E12000002	E12000004	North West	East Midlands	259
E12000008	E12000004	South East	East Midlands	315
E12000009	E12000004	South West	East Midlands	72
W92000004	E12000004	Wales	East Midlands	30
E12000005	E12000004	West Midlands	East Midlands	841
E12000003	E12000004	Yorkshire and The Humber	East Midlands	751

Of these, we can get the total flows emitted from every origin \(O_i\) and equivalently the arriving at \(D_j\):

sum_dest <-  aggregate(t_wu02Ew_od$OBS_VALUE,
                by = list(t_wu02Ew_od$USUAL_RESIDENCE_CODE),
                FUN = sum)
names(sum_dest) <- c('Destination','Dj')

sum_orig <-  aggregate(t_wu02Ew_od$OBS_VALUE,
                by = list(t_wu02Ew_od$PLACE_OF_WORK_CODE),
                FUN = sum)
names(sum_orig) <- c('Origin','Oi')

Bring all the above together we get:

t_wu02Ew_od <- merge(t_wu02Ew_od, dist_pl, by.x=c("USUAL_RESIDENCE_CODE", "PLACE_OF_WORK_CODE"), by.y=c("Origin", "Destination"))

t_wu02Ew_od <- rbind(t_wu02Ew_od, sum_dest[sum_dest$Origin %in% t_wu02Ew_od$PLACE_OF_WORK_CODE])
t_wu02Ew_od <- merge(t_wu02Ew_od, sum_dest, by.x="PLACE_OF_WORK_CODE", by.y= "Destination")
t_wu02Ew_od <- merge(t_wu02Ew_od, sum_orig, by.x="USUAL_RESIDENCE_CODE", by.y= "Origin")

kable(head(t_wu02Ew_od))

USUAL_RESIDENCE_CODE	PLACE_OF_WORK_CODE	USUAL_RESIDENCE_NAME	PLACE_OF_WORK_NAME	OBS_VALUE	Distance_m2	Dj	Oi
E12000001	E12000001	North East	North East	18512	0.0	19158	19262
E12000001	E12000004	North East	East Midlands	44	243896.0	46090	19262
E12000001	E12000008	North East	South East	47	425126.6	107425	19262
E12000001	E12000007	North East	London	79	409329.8	74498	19262
E12000001	E12000002	North East	North West	75	121184.5	67584	19262
E12000001	E12000006	North East	East	34	348467.2	71042	19262
Run the model:

#do not consider within flows
t_wu02Ew_od <- t_wu02Ew_od[t_wu02Ew_od$Distance_m2 != 0,]

#note: I have not applied a function to penalise for increasing distance. This is just for simplicity of the example
unconstrained_gm <- glm(OBS_VALUE ~ log(Oi)+log(Dj)+log(Distance_m2), na.action = na.exclude, family = poisson(link = "log"), data = t_wu02Ew_od)
#let's have a look at it's summary...
summary(unconstrained_gm)

## 
## Call:
## glm(formula = OBS_VALUE ~ log(Oi) + log(Dj) + log(Distance_m2), 
##     family = poisson(link = "log"), data = t_wu02Ew_od, na.action = na.exclude)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -60.133  -11.135   -5.893    1.386  100.783  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      11.76801    0.32101   36.66   <2e-16 ***
## log(Oi)           0.61235    0.01643   37.26   <2e-16 ***
## log(Dj)           0.56392    0.01572   35.87   <2e-16 ***
## log(Distance_m2) -1.54905    0.00801 -193.40   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 94447  on 71  degrees of freedom
## Residual deviance: 23912  on 68  degrees of freedom
## AIC: 24432
## 
## Number of Fisher Scoring iterations: 5

rsq(unconstrained_gm)

## [1] 0.5649217

Given that the data is in matrix form, there are ways of taking into account the row and column totals to calibrate the model. At a higher level, these are a ratio (or percent) for each cell, which takes into account the differences between the observed row, then column totals compared to what we would get if we only considered the distance (or more generally cost) associated with that particular type of flow. Note, I say flow, not cell, because we can assume that the cost of moving from some groups of places might be associated with the same cost and would therefore use the same calibrating ratio.

Wolwe, Buragard and Breisslein go into this in more detail, whilst discussing fairly new R package for trade gravity modelling.

Further reading: Gravity models of Trade in R
Unconstrained Model of Migration